The invention relates to systems for solving finite element method (FEM) equations.
The equilibrium governing equation in the finite element method (FEM) may be written as follows: EQU KU=R, (1)
where K is a stiffness matrix, U is a displacement vector specifying the unknown nodal displacements, and R is a load vector. The load vector includes the effect of element body forces, the element surface forces and the element initial stresses and concentrated loads on the assemblage of elements or nodes. The stiffness matrix describes the structural properties and inter-element interactions within the body of elements being characterized. In the FEM, the calculation of K is accomplished by integrating material properties over a set of nodes. Descriptions of the commonly used techniques for calculating K are found in many standard texts on finite element analysis. (See, for example, J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice-Hall, 1982 and L. J. Segerlind, Applied Finite Element Analysis, John Wiley and Sons, 1984.)
There are many ways to solve Eq. 1, however, because of the typically large size of the K matrix direct integration methods are preferred. For example, a first-order iterative solution technique is commonly used: ##EQU1## where .DELTA.t is the integration time step and c is the system damping factor. This solution is equivalent to direct integration of EQU CU+KU=R, (3)
where C is a diagonal matrix whose entries are the damping factor c.
A much more efficient approach to solving Eq. 1 is to employ a second-order solution technique such as the Central Difference Method (CDM) for dynamic analysis. Such a second order solution technique is equivalent to direct integration of: EQU MU+CU+KU=R (4)
for some n.times.n matrices M and c.
In the FEM, the matrix M is interpreted as a mass matrix describing the distribution of mass within the body, and c is interpreted as a damping matrix describing the body's internal dissipation of energy. The process of direct integration is interpreted as stepping forward by fixed time increments to simulate the reaction of the body to a load R. The equilibrium nodal displacements U describe the shape of the body after it has come to rest, and are the solution to the initial loading problem given by Eq. 1.
Referring to Eq. 4, to obtain an equilibrium solution U, one integrates the equation using an iterative numerical procedure (such as CDM) at a cost of roughly 3nm.sub.k operations per time step, where n is the order of the stiffness matrix and m.sub.k is its half bandwidth. (See Bathe, Appendix A.2.2 and Segerlind for complete discussions on the bandwidth of a stiffness matrix.) Thus, there is a need for a method which transforms the Eq. 4 into a form which leads to a less costly solution. Since the number of operations is proportional to the half bandwidth mk of the stiffness matrix, a reduction in m.sub.k will greatly reduce the cost of step-by-step solution.
To accomplish this goal a transformation of the nodal point displacements U can be used: EQU U=PU, (5)
where P is a square transformation matrix and U is a vector of generalized displacements. Substituting Eq. 5 into Eq. 4 and premulti by P.sup.T yields: EQU MU+CU+KU=R (6)
where EQU M=P.sup.T MP; C=P.sup.T CP; K=P.sup.T KP; R=P.sup.T R.
With this transformation of basis set a new system of stiffness, mass and damping matrices can be obtained which has a smaller bandwidth then the original system.
The optimal transformation matrix is derived from the free vibration modes of the equilibrium equation. This can be seen by examining the undamped case. Under the assumption of no damping, the governing equation reduces to: EQU MU+KU=R (7)
From this, an eigenvalue problem can be derived which will determine an optimal transformation basis set .phi.: EQU K.phi.=.omega..sup.2 .phi.M (8)
The eigenvalue problem in Eq. 8 yields n eigensolutions: EQU (.omega..sub.1.sup.2,.phi..sub.1),(.omega..sub.2.sup.2,.phi..sub.2), . . . ,(.omega..sub.n.sup.2,.phi..sub.n),
where all the eigenvectors are M-orthonormalized. Hence, ##EQU2## and EQU 0.ltoreq..omega..sub.1.sup.2 .ltoreq..omega..sub.2.sup.2 .ltoreq..omega..sub.3.sup.2 .ltoreq.. . . .ltoreq..omega..sub.n.sup.2( 10)
The vector .phi..sub.i is the i.sup.th mode shape vector and .omega..sub.i is the corresponding frequency of vibration. Now I define a transformation matrix .PHI., which has for its columns the eigenvectors .phi..sub.i, and a diagonal matrix .OMEGA..sup.2, with the eigenvalues .omega..sub.i.sup.2 on its diagonal: ##EQU3## Equation 8 can now be written as: EQU K.PHI.=.OMEGA..sup.2 .PHI.M, (12)
and since the eigenvectors are M-orthonormal: EQU .PHI..sup.T K.PHI.=.OMEGA..sup.2 ( 13) EQU .PHI..sup.T M.PHI.=I. (14)
From the above formulations, it becomes apparent is the optimal transformation matrix P for systems in which damping effects are negligible, since for this transformation the equilibrium equation is reduced to EQU U+.OMEGA..sup.2 U=.PHI..sup.T R(t) (15)
or, equivalently, to n independent and individual equations of the form EQU U.sub.i (t)+.omega..sub.i.sup.2 U.sub.i (t)=r.sub.i (t) 916)
where r.sub.i (t)=.PHI..sub.i.sup.T R(t) for i=1,2,3, . . . , n.
Applying this transformation to the equilibrium governing equation (Eq. 1) yields the following results: EQU K.PHI.U=R. (17)
Multiplying both sides of Eq. 17 by .PHI..sup.T : EQU .PHI..sup.T K.PHI.U=.PHI..sup.T R, (18)
which, using Eq. 13, becomes: EQU .OMEGA..sup.2 U=.PHI..sup.T R, (19)
or, equivalently: EQU U=.OMEGA..sup.-2 .PHI..sup.T R. (20)
To get back to the original nodal point displacement reference system, the transform .PHI. is applied to both sides of Eq. 20: EQU .PHI.U=U=.PHI..OMEGA..sup.-2 .PHI..sup.T R=W R, (21)
where W=.PHI. .OMEGA..sup.-2 .PHI..sup.T.
In other words, using transformation .PHI. yields a closed form solution that can be computed by multiplying R by W. Since W may be precomputed, the number of operations required at the time of computing U is O(n.sup.2), a significant improvement over the O(n.sup.3) operations required for the previously described integration methods. Nevertheless, this technique still requires a very substantial amount of precomputation. More seriously, however, whenever discontinuities are introduced into the system the stiffness matrix w must be recomputed. Even if perturbation methods are employed to recompute W, the computational costs can be significant. It is desirable, therefore, to find a simpler, more efficient method of calculating W, of performing the multiplication of W and R, and of compensating for the introduction of discontinuities.